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SUMMARY 


New flux-corrected transport algorithms are described for solving gener- 
alized continuity equations. These techniques were developed by requiring 
that the finite-difference formulae used ensure positivity for an initially 
positive convected quantity. Thus FCT is particularly valuable for fluid-like 
problems with strong gradients or shocks. Repeated application of the same 
subroutine to mass^ momentum^ and energy conservation equations gives a simple 
solution of the coupled time-dependent equations of ideal compressible fluid 
dynamics without introducing an artificial viscosity. FCT algorithms span 
Eulerian, sliding-rezone, and Lagrangian finite-difference grids in several 
coordinate systems. The latest FCT techniques are fully vectorized for 
parallel/pipeline processing and are being used on the Texas Instruments ASC 
at NRL. 


INTRODUCTION 


This paper reviews the Flux-Correct Transport (FCT) 
been developed to solve the continuity equation 


techniques which have 
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which get reflected in some of the numerical solution techniques: 
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The convective term V*^p displayed explicity in Eq. (lb) gives (la-c) its 
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intrinsically hyperbolic form and causes really severe problems numerically. 
The compression term - is sometimes absent^ as in the Liouville 

Equation^ in which case it is often called the convection or the advection 
equation. 

Continuity equations underlay compressible and incompressible fluid 
dynamics^ hydrod 3 mamicS; plasma physics (Vlasov Equation and MHD moment 
equations) and even quantum mechanics. They appear in most descriptions of 
dynamic physical systems simply because they express two of the more general 
principles in physics^ conservation and causality. Continuity equations also 
display the positivity property: a quantity being transported will never turn 
negative anywhere in a reasonable flow field if that quantity was everwhere 
positive to start with. This positivity property expresses in a continuum 
way the intrinsic corpuscular nature of matter. Thus matter cannot be re- 
moved from a region which is devoid of matter to begin with. This duality^ 
in which matter obeys both particle and fluid-like equations on microscopic 
and macroscopic scales respectively^ has its ramifications for numerical 
solution techniques as well. The first three of the possible solution tech- 
niques listed below take advantage of the underlying discrete basis of the 
continuity equation while the last three aim more directly at solving the 
partial differential equation itself. 

TABLE 1 - Possible Solution Techniques 

• Quasiparticle Methods 

1. Collisionless particles^ stars and plasmas 

2. Collisional particles for fluids 

• Charac te r i s t ic Me thod s 

• Lagrangian Finite Difference Methods 

• Eulerian Finite Difference Methods 

1. Explicit vs implicit 

2. Order vs accuracy 

• The Finite Element Method 

• The Spectral Method 

Because of their speed and simplicity^ finite-difference solutions of 
the continuity equation must always be considered carefully as the most likely 
of the six candidate methods (ref l). Quoting three conclusions in a similar 
context from another source (ref 2)*. 

• The spectral method has no stability problems but is much more 
complicated and slower than generalized difference methods. 

• It is doubtful whether the finite-element method^ based on piece- 
wise polynomials^ can compete with the above methods. 
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• If difference methods are used^ they should be at least fourth- 
order accurate. 

While I agree basically with these remarks^ they do not really encompass the 
quasiparticle schemes nor do they adequately reflect a very important piece 
of personal experience. Whenever a theory is to be tested^ or an algorithm^ 
or a new mathematical technique^ attention rather quickly turns to the crucial 
yet simple conservation equations of ideal compressible flow. 

Finite-difference methods have solved the transient Rankine-Hugoniot 
shock problem adequately. For that matter^ various flavors of quasiparticle 
methods can do the same vital problem creditably if enough particles are used. 
I do not know of any calculation^ even in one spatial dimension, using either 
a finite-element or a spectral method which has correctly solved for an ideal 
gas compressible shock. Until such calculations become common place and com- 
putationally attractive, finite-difference methods would seem to have the 
inside track. 


IMPROVING FINITE-DIFFERENCE TECHNIQUES 

There is undoubtedly some merit in trying to boost the performance of 
quasiparticle methods on the one hand and the basis -function expansion methods 
on the other toward the performance obtained from finite differences. But 
it is a low risk-high return investment to patch up the obvious failings of 
the front runners, finite differences. In the case of Lagrangian finite- 
difference methods, the major outstanding problems arise from secularly un- 
attractive distortions of the grid which wreck calculations of interesting 
flows quite quickly (refs. In the case of Eulerian methods, the major 

outstanciing weakness in a hugh class of problems of real interest is the need 
for a large artificial damping (numerical diffusion) to fill in what would 
otherwise be pits of ^negative density” in the calculated profiles. Since the 
”Eulerian" positivity problem is encountered even in Lagrangian calculations 
for many situations, it demands the greater share of attention. 

The Flux-Corrected Transport techniques (FCT) which are the subject of 
this paper have been designed carefully to satisfy the following six require- 
ments of an "ideal” algorithm for solving the continuity equation (ref, 6-9)* 
An ideal algorithm should: 

1 . Be linearly stable for all cases of interest, 

2 . Mirror conservation properties of the physics, 

5- Ensure the positivity property when appropriate, 

4. Be reasonably accurate, 

5. Be computationally efficient, and 

6. Be independent of specific properties of one application* 
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FCT algorithms arise naturally (ref 6) as a result of trying to satisfy 
requirement 3* Most work in the past has centered on trying to increase the 
mathematical order of accuracy of a scheme while ignoring the physical 
positivity property which the fluids display prominently. 


Consider the rather general three-point approximation to Eq. (la) 



where 




is the density at mesh point j. 


Equation (2) is in finite-difference conservative form with whole indices 
representing cell centers and half indices indicating cell interfaces. The 
additional numerical diffusion terms with diffusion coefficients '^•-(- 1/2 


have to be added to ensure positivity. The stability of Eq. (2) is ensured^ 
at least roughly^ when 



^j+l/2 



-2 

j+1/2 


(5a) 


The upper limit arises from the explicit diffusion time-step condition while 
the lower limit is the Lax-Wendroff damping. Unfortunately positivity is 
only ensured linearly when 




> 1 
j+l/2 2 


®j+l/2 


the firs t-order> ups tream-centered scheme result. 


(5b) 


We appear to be caught between a rock and a hard place here but the 
escape route is signaled in the preceding sentence by the word "linearly". 

By relaxing the linearity implied by Eq. (2) and letting the diffusion 
coefficients be nonlinear functionals of the flow velocities <®j-|-i/2l ^ 

we can hope to reduce the integrated dissipation below the ramer ghastly 
limit (5b) and yet retain sufficient dissipation near steep gradients to 
ensure positivity. A literature is beginning to form about these "monotonic" 
difference schemes (refs. 6-9^ 10, 11 ) since the dilemma of accuracy versus 
positivity in Eulerian difference schemes can be resolved in no other way. 
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FLUX-CORRECTED TRANSPORT ALGORITHMS 


The firsts and so far the most developed and used^ of the monotonic 
schemes is Flux-Corrected Transport. The calculation in Fig. 1 was performed 
by the first FCT Algorithm SHASTA (ref. 7) and had an error about four or 
five times smaller than the simple linear methods also shown. The damping 
was second order as were the relative phase errors. Figure 1 shows a 
comparison of four common difference schemes solving the standard square wave 
problem. The effects of excess numerical damping in the donor-cell treatment 
(upstream-centered first-order), and of excess dispersion in the leap frog 
and Lax-Wendroff treatments are clearly visible. Dispersion manifests itself 
as a trail or projection of oscillations in the computed solution near dis- 
continuities and sharp gradients of the "correct** solution. 

The basic FCT technique shown in Figure 1 was quickly generalized to 
cylindrical and spherical systems, to Lagrangian as well as fixed Eulerian 
grids, and was applied to a number of one-, two-, and three-dimensional 
problems. More recent work has been devoted toward extending the basic non- 
linear flux- correct! on techniques to convection algorithms other than SHASTA 
and toward discovering an **optiraum** FCT algorithm. 

Since the latest FCT algorithms have eliminated roughly 95^ of the 
removab le error and the removable error that remains is barely half of the 
irreducible error, it is natural to have turned next toward optimization 
in speed, flexibility, and generality (ref. 12). In Flux-Corrected Transport 
algorithms, the basic convective transport algorithm is augmented with a 
strong enough linear diffusion to ensure positivity at the expense of excess 
smoothing. Since the amount of diffusion which has been added is known, FCT 
then performs a conservative antidiffusion step to remove the diffusion in 
excess of the stability limit. However, the antidif fusive fluxes are 
effectively multiplied by a coefficient which ranges from zero to unity to 
preserve monotonicity. The criterion for choosing the reduction factors of 
the antidif fusive fluxes is that the antidiffused solution will exhibit no new 
maxima or minima where the diffused solution had none. 


Although the limit (5b) represents the minimum amount of diffusion 
needed for stability, FCT algorithms generally use a larger zero-order 
diffusion because it has been found that the correct choice of the 

within the mono tonic and stable range will reduce convective phase 




errors xrom second to fourth order as suggested by Kreiss. Since the anti- 
diffusion can also be chosen correspondingly larger, no real price is exacted 
for this improvement in phase properties. The most recent efforts (ref. 12) 
have taken advantage of this fact to generate minimum- operation-count FCT 
algorithms for Cartesian, cylindrical, and spherical coordinate systems with 
stationary (Eulerian) and moveable (Lagrangian) grid systems. Because these 
algorithms, and in particular the nonlinear flux-correction formula, were 
very carefully designed, they are fully **vectorizeable** for pipeline and 
parallel processing and have been implemented in all generality in Fortran on 
the Texas Instruments Advanced Scientific Computer at the Naval Research 
Laboratory. The execution time per continuity equation per grid point is 


1295 


roughly I .5 |j,sec. A complete 2D calculation on a system of 200 grid points 
X200 grid points requires roughly 2 to 2.5 seconds per timestep depending on 
extra physics and boundary conditions incorporated in the problem. 

Figure 2 shows a ID calculation performed by the code FASTID on the ASC. 

The problem chosen is the "Lapidus” problem (ref. I 3 ) in cylindrical coordi- 
nates with y = 1.4. The diaphragm is originally at r - 1.0 in Fi& 2 and 
bursts at t - 0.0 sec. The density solution is shown at 0.6 sec in Fig. 2 
just after the shock has reached the axis and rebounded. Three different 
resolution calculations are overlapped to show the convergence. Even the 
calculation with 50 cells is at least as accurate as the original Payne and 
Lapidus solutions with 200 cells. The width of the contact discontinuity 
is about 5*5 cells while the shock is smeared over only I.5 cells without 
noticeable overshooting or undershooting. The calculations are performed without 
any added artificial viscosity. The monotonicity control provided by FCT 
on each of the continuity equations separately is adequate to ensure stability 
and accuracy as shown. 

Figures Ja^ b show the evolution of Rayleigh-Taylor instability in the 
implosion of a laser pellet shell 50 microns thick. The calculation was 
performed using the FAST2D code on the ASC and the thermal conductivity 
was set to zero to show the full nonlinear evolution of the instability. An 
Eulerian "sliding-rezone" grid was used with 200 X 200 grid points. Strong 
deterioration of the shell is apparent by 3.85 nsec but breakthrough has not 
yet occurred. The jetting of material off the backside is severe enough by 
this time that grid distortion would obviate the usual Lagrangian solution 
techniques. Figure shows the linear phase of the instability and shows 
the onset of the nonlinear regime. 
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Figure 1.- Comparison of four difference schemes solving 
the square wave problem. The merits of Flux-Corrected 
Transport relative to the other methods are clear. 
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Figure 2.- A cylindrical diaphragm problem using FASTID. Three different resolutions 
are shown to demonstrate convergence of the solution. FCT generally requires shock 
widths of about 1.5 cells. 
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of the Rayleigh-Taylor instability. The shell is 30 microns thick with 



